Measurement of the reduced scattering 
coefficient of turbid media using single 
fiber reflectance spectroscopy: fiber 
diameter and phase function dependence 

S. C. Kanick, 1 2 * U. A. Gamm, 1 2 M. Schouten, 1 
H. J. C. M. Sterenborg, 1 D. J. Robinson, 1 and A. Amelink 1 

1 Center for Optical Diagnostics and Therapy, Department of Radiation Oncology, Erasmus 
Medical Center, PO Box 2040, 3000 CA Rotterdam, The Netherlands 
^Denotes that authors contributed equally. 
*s. kanick@erasmusmc.nl 

Abstract: This paper presents a relationship between the intensity 
collected by a single fiber reflectance device {Rsf) and the fiber diameter 
(dfib) and the reduced scattering coefficient (/i') and phase function (p(0)) 
of a turbid medium. Monte Carlo simulations are used to identify and 
model a relationship between R$f and dimensionless scattering (jj.[dfjb). 
For n' s dfjt, > 10 we find that R$f is insensitive to p(0). A solid optical 
phantom is constructed with fl' s s» 220 mm~ 1 and is used to convert Rsf of 
any turbid medium to an absolute scale. This calibrated technique provides 
accurate estimates of \i[ over a wide range ([0.05 — 8] mm ) for a range of 
d fib ([0.2-1] mm). 
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1. Introduction 

Reflectance spectroscopy is widely used for noninvasive measurement of the tissue absorp- 
tion and scattering coefficients, which can provide diagnostic information about tissue function 
and structure, respectively [1-9]. Scattering in tissue originates from spatial heterogeneities of 
the optical refractive index that occur on size scales ranging from a few nanometers to a few 
millimetres [9]. Scattering properties of tissue include the scattering coefficient (/i s ), which de- 
scribes the mean free path between scattering events, and the scattering phase function (p(9)), 
which describes the angular distribution of scattering events. These properties are often de- 
scribed in terms of the scattering anisotropy g, which is the expectation value for the cosine 
of the scattering angle g =< cos(6) >, and the reduced scattering coefficient jj,' s = jU s (l — g), 
which is a description of the combined effect of scattering coefficient and average scattering 
angle. Reflectance spectroscopy measurements with large source-detector separations (« 10 
mean free paths [10]) can be modeled using diffusion theory, which utilizes /i( to character- 
ize the effect of scattering on reflectance. Diffuse measurements return information about bulk 
tissue properties, which is a limitation because many clinical diagnostic applications require 
more localized measurements [11-13]. Quantitative extraction of the absorption and reduced 
scattering coefficients from reflectance measurements with small source-detector separations is 
complicated because light transport in this non-diffuse regime, where collected photons may 
undergo few scattering events, is no longer dependent exclusively on jj.' s , but is also sensitive to 
large angle back scattering events, a characteristic of the scattering phase function. 

The influence of phase function on reflectance intensity measured close to the source has 
been described in the literature. For example, Mourant et al. [14] reported significant differ- 
ences in the Monte Carlo simulated reflectance intensities for a small source-detector probe for 
different phase functions (Henyey-Greenstein vs. Mie) but with the same anisotropy. Kienle et 
al. [15] also reported large deviations in the calculation of jj. a and fx' s after simulating the re- 
flectances for a source-detector separation of < 1 mm for multiple phase functions by applying a 
standard solution for the diffusion equation [10] to extract the optical properties. Bevilacqua and 
Depeursinge [16, 17] systematically investigated the influence of higher order moments of the 
phase functions for several source-detector separations on the reflectance intensities. They in- 
troduced a phase function dependent coefficient y which, together with /i„ and [i' s and refractive 
index n, would give a complete description of the reflectance for systems with measurements 
from multiple short source detector separations. Hull and Foster later developed an analyti- 
cal description of light transport termed the P3 approximation, which models both diffuse and 
non-diffuse light scattering by including the first three moments of the phase function [18]; this 
approach has been shown to be valid for source-detector separations as low as f=a 0.5 mm [13], 
but not directly applicable for devices with overlapping source and detector locations. Other 
empirical investigations of the relationship between reflectance intensity with small source- 
detector separations and /ij [20-22] did not fully characterize the effect of phase function on 
the collected signal. 
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Recently, our group reported on the use of experimentally validated Monte Carlo (MC) sim- 
ulations to investigate the effect that optical properties have on the propagation of photons that 
are collected during a single fiber reflectance (SFR) measurement [23]. SFR spectroscopy uses 
a single fiber for the delivery and detection of light; a schematic of the device setup is shown 
in Fig. 1. The advantages of an SFR device are the simplicity of the setup machinery and the 
small dimensions of the fiber-optic probes, which allows e.g. measurement through the lumen 
of FNA-needles [24, 25]. In the past, single fiber geometries have been used for the investi- 
gation of fluorescence in tissue and tissue simulation phantoms [26-29] and for particle size 
analysis [30, 31]. The relationship between single fiber reflectance intensity and optical prop- 
erties was studied by Moffit and Prahl [19]; however, their analysis did not yield quantitative 
estimates of fl' s or characterize phase-function dependence of the signal. Further studies have 
focused on technical aspects of single fiber reflectance, such as the need to remove internal 
specular reflections from the collected signal [32] and the dependence of the collection effi- 
ciency of a single fiber on the optical properties of the sampled medium [33]. Our previous 
work has focused on utilizing SFR spectroscopy to quantitate the absorption coefficient (/i„), 
and the corresponding chromophore concentrations, in turbid media. These studies introduced 
an empirical function that describes the dependence of the effective photon path length [23, 34] 
on the optical properties of the sampled medium. A limitation of our method is that the path 
length relation requires knowledge of jj.' s , a parameter that is currently not quantitatively ex- 
tracted from the measured spectra. The current spectral analysis algorithm uses estimated values 
for the wavelength-dependent jj.' s and the uncertainty introduced has been documented [24] . A 
quantitative measurement of jj.' s would improve the accuracy of our chromophore concentration 
estimates. Moreover, quantitative measurement of fl' s using an SFR device is inherently useful, 
because it may provide information about the tissue microarchitecture [9]; information that is 
complementary to the physiological information extracted from jj. a , and may have diagnostic 
value. Additionally, quantitative knowledge of jj.' s and extracted from SFR measurements 
may be used to correct for the effect of optical properties on single fiber fluorescence measure- 
ments [26-29]. 




Bifurcated Fiber 



Single Fiber 



Fig. 1. Schematic of single fiber reflectance probe machinery. 



There are two challenges associated with quantitative extraction of /i' from an SFR measure- 
ment. First, a calibration procedure must be designed that allows absolute measurement of the 
reflectance (Rsf) of a turbid medium from an SFR measurement. Second, Rsf must be related 
to /J.' s . The first challenge will be addressed by designing a calibration phantom that returns a 
known percentage of incident light during measurement. This measurement allows the scaling 
of reflectance intensity measured on a turbid medium to a percentage of the absolute amount 
of incident measurement light. Although the concept of using calibration phantoms to scale the 
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reflectance is not novel, for an SFR device (or any other small source-detector separation de- 
vice) the scattering phase function of the calibration phantom is a critical factor that needs to be 
accounted for in the design of the phantom; the present study considers phase function effects 
in detail. The second challenge is addressed by investigating the relationship between R$p and 
the scattering properties of the medium, and the influence of scattering phase function using an 
extensive set of Monte Carlo simulations. The simulations are experimentally validated using 
highly scattering ( fi' s w 20 — 200 mm~ 1 ) solid phantoms consisting of titanium dioxide (T1O2) 
suspended in silicone, as well as low and medium scattering (0.05 < fl' s < 8 mm -1 ) liquid phan- 
toms consisting of a range of dilutions of Intralipid. A semi-empirical model is introduced that 
describes the relation between R$f and dimensionless scattering (jj.' s dfit) for a wide range of 
phase functions. 

Our methodology to experimentally extract /i' of a turbid medium from SFR spectra, us- 
ing a calibration phantom and a mathematical model of R$f = f(n' s dfn,), is validated using 
SFR measurements on Intralipid phantoms. We show that with knowledge of the wavelength- 
dependent phase function, the model function can be used to accurately extract from SFR 
measurements for a wide range of reduced scattering coefficients (jl' s = [0.05 — 8] mm ) and 
fiber diameters (d fib = [0.2,0.4,0.6,0.8,1] mm). 

2. Materials and methods 

2.1. Monte Carlo simulation 

The Monte Carlo (MC) code used to simulate single fiber reflectance spectroscopy in this study 
has been described in detail previously [23]. In short, our MC code is based on the MCML 
code [35] that stochastically simulates photon propagation within a turbid medium. During pho- 
ton propagation, each photon step size was selected from an exponential distribution weighted 
by the scattering coefficient, and each scattering angle was selected from a user-specified phase 
function. Reflection and refraction due to the index of refraction mismatch at the medium/fiber 
and medium/air interface were calculated using the Fresnel equations and Snells law. The in- 
dex of refraction («) of the medium and fiber were specified at 1.37 and 1.5, respectively. The 
numerical aperture (NA) of the fiber was set as 0.22. Photons were initialized by selecting a 
location from a uniform distribution on the single fiber face in contact with the turbid medium 
(z = 0), and the launch direction was selected from a uniform distribution of angles within the 
fiber cone of acceptance, where the acceptance angle was given as Q a = asin ( — ) . For 

\ n medium / 

photons propagating within the turbid medium that cross the medium interfacial boundary at 
z = 0, it was checked if they hit the fiber face; those in contact with the fiber face and traveling 
at an angle within the fiber cone of acceptance were collected, the rest were terminated. Pho- 
tons propagating within the medium far from the fiber face do not contribute to the collected 
reflectance intensity and were terminated at a hemispherical limit from the fiber face of 10-^r-; 
a limit that was confirmed to not influence model outputs for the range of optical properties 
investigated in this study. For each simulation, the MC code outputted the reflectance inten- 
sity as the ratio of the total number of photons collected (TPC) to the total number of photons 
launched (TPL) during the MC simulation, calculated as: 

R MC _ TP C q\ 
SF TpL 

MC simulations were run over a broad range of reduced scattering coefficient values (jiij = 
[0.1,0.3,1,3,10,30,100,300] mm" 1 and 5 different fiber diameters [0.2,0.4,0.6,0.8,1.0] mm. 
To investigate the dependence of the collected signal on phase function, sets of simulations were 
run that investigated variations in both phase function and anisotropy values. Figure 2 shows the 
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the angular scattering probability of selected phase functions. The Henyey-Greenstein-phase 
function (HG), commonly used to describe photon scattering in tissue, was simulated for g- 
values of [0.5,0.7,0.8,0.9]. The Modified Henyey-Greenstein phase function (MHG), which 
contains a pronounced backscattering feature, was simulated with an effective anisotropy value 
of g = 0.9 as in [21,23]. In order to compare simulations with experimental measurements in 
Intralipid 20%, simulations were performed with the wavelength-dependent phase functions of 
Intralipid as given by Michaels et al. [36]. Simulations were run for wavelengths from 400 to 
900 nm in steps of 50 nm; corresponding anisotropy values for these phase functions were in 
the range of [0.47 — 0.82]. At least 2 million photons were launched for each MC simulation. 

icr 2 



1CT 3 



i "itr 4 

a. 

10" 6 



0 30 60 90,_, , 120 150 180 

Fig. 2. Angular distribution of scattering events for selected scattering phase functions. 




2.2. Semi-empirical model of the single fiber reflectance intensity 

Inspection of simulated reflectance data over the range of simulated scattering properties 
yielded an observable relationship between R¥£ and the product fl' s dfn,\ this product is re- 
ferred to throughout this paper as dimensionless scattering. A semi-empirical model is utilize 
to describe this relationship, given by the following set of Equations, 



R 



Model 
SF 



TPC 
TPL 
TPH 
TPL 
TPC 
TPH 



P: 



(2) 
(3) 
(4) 



The collected reflectance intensity (Rfp del ) is calculated as the ratio of TPC to TPL during the 
MC simulation, as given in in Eq. (2). This ratio can be described by the the product of the 
fraction of photons remitted from the medium that contact the fiber face (O), and the collection 
efficiency (t] c ) of the fiber. Equation (3) presents O as the ratio of photons that are remitted 
from the medium and contact the fiber face (TPH) at all angles vs. the TPL. The right hand 
side of Eq. (3) introduces a mathematical expression that captures a saturating relationship [37] 
between O and dimensionless scattering that was observed in the MC outputs. The collection 
efficiency (rj c ) is given in Eq. (4) as the ratio of the photons of TPCITPH, which represents 
the fraction of all photons remitted from the tissue that contact the fiber face within the cone 
of acceptance of the fiber. The concept of r\ c has been described in detail previously by Bargo 
et al. [33]; in brief, t] c is dependent on optical properties but for high scattering coefficients 
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it collapses to an approximate lower limit of rj/,,,,,-, ~ sin 2 (0„), where the acceptance angle is 
dependent on fiber NA and the refractive index of the medium, as 0„ = asin(NA / n me <a um ) . 
The observed dependence of both O and rj c on dimensionless scattering is presented in Section 
3. 1 .2 and described in Section 4. 

The parameters [p i , p2 , P3] were fitted by minimizing the weighted residual error between the 
simulated (Rff) and model-estimated (R^p del ) reflectance intensities, with each point weighted 
by the inverse of the simulated data point. It is important to note that neither O nor 7] c are 
observable from our experiments, and while they were estimated by the MC model, the fitted 
parameter values [pi,p2,p3] were only fitted to reflectance intensity data using Eq. (2); i.e. 
we did not separately fit Eqs. (3) and (4) which represent unobservable quantities. Parameters 
were estimated using a Levenberg-Marquardt algorithm [38] written into Matlab code (version 
R2009a, MathWorks). Confidence intervals of the parameter estimates were calculated from 
the square root of the diagonal of the covariance matrix, as described in detail elsewhere [39]. 

2. 3. Experimental setup for single fiber reflectance spectroscopy 

Figure 1 shows a schematic of the single fiber reflectance spectroscopy device setup, which 
consists of a halogen light source (HL-2000-FHSA, Ocean Optics, Duiven, Netherlands), a 
spectrometer (SD2000, Ocean Optics, Duiven, Netherlands) and a solid core single fiber probe. 
The single fiber is connected to the light source and the spectrometer via a bifurcated fiber. Dur- 
ing measurements, light travels from the light source through the fiber and exit into the probe 
tip that is in contact with the sampled medium. Photons that backscatter to the fiber face at an 
incident angle that is within the cone of acceptance enter the fiber core and travel to the spec- 
trometer. To remove specular reflections within the collected reflectance intensity due to index 
of refraction mismatch at the probe/medium interface tip, the probe tip is polished at an angle 
of 15 degrees [34]. Variations in the output of the lamp, transmission characteristics of the fiber, 
and sensitivity of the spectrometer, as well as remaining specular reflectance and other internal 
reflections are being accounted for by performing a calibration procedure that we have de- 
scribed previously [34]. The calibration includes measurements of white and black spectralon 
standards (Labsphere SRS-99 and SRS-02, spectrally flat, with 99% and 2% reflectance, re- 
spectively) and a measurement of water in a black container. The single fiber reflectance signal 
can be calculated as 

* I water 
_I white ' 



K SF 



(5) 



Here I wa ter, Iwhite, an d hlack are the measured intensities of water, and black and white spectralon 
respectively and / stands for the measured intensity of the sample. Because the signals I w hite> an d 
hlack depend on the distance of the fiber probe to the spectralon during calibration measurement, 
jfmeas-rel ca i cu i a t ec [ m gq (5) j s relative to the spectralon distance. 

2.4. Solid phantom preparation 

Solid phantoms were prepared according to the recipe reported by de Bruin et al. [40]. Silicone 
[Sylgard 184 Silicone Elastomer Kit (Dow Corning Europe SA, Seneffe, Belgium)] served as 
a stable matrix and Titanium dioxide (Ti02) (Sigma Aldrich, Zwijndrecht, Netherlands) was 
used as a scatterer. All solid phantoms were prepared with a weight of 20 g, including different 
amounts of Ti02 to create a series of different highly scattering phantoms (Ti02 weight percent - 
age= [0,2.5,5, 10,20,30] %). The silicone kit consisted of the silicone compound and a curing 
agent which was always mixed in a ratio of 9:1. To prepare the phantoms the silicone com- 
pound and Ti02 were weighed, mixed with a Dremel 300 (Dremel, Leinfelden-Echterdingen, 
Germany) and placed into an ultrasound bath for 10 min, to assure that the Ti02 particles were 
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homogeneously distributed within the silicone. Then the curing agent was added and the com- 
position was mixed again until the agent was completely mixed into the compound. The mass 
was poured into the mold and placed in a vacuum desiccator (Duran Group GmbH, Mainz, 
Germany) and held under vacuum until remaining air bubbles were removed from the phantom 
(1—3 hours). The phantoms were cured after approximately 12 hours at room temperature. 

During the single fiber reflectance measurements the single fiber probe tip was brought into 
contact with the phantom surface, with no gap between fiber face and phantom surface. Three 
spectra were recorded at 3 different locations on the phantom surface and averaged. Each phan- 
tom was measured with the 5 different fiber probes with diameters of [0.2,0.4,0.6,0.8,1.0] 
mm. To test the homogeneity of the phantoms and the reproducibility of the preparation, 6 
solid phantoms with a weight percentage of 30 % Ti02 were prepared and single fiber spectra 
were measured at 10 different locations on each phantom surface. The maximum relative devi- 
ation between measurements on two phantoms was 0.66 %, while the intra-phantom variability 
was even lower. 

2.5. Liquid phantom preparation 

The liquid phantoms were prepared by mixing 0.9 % NaCl (Baxter, Utrecht, Netherlands) 
with different amounts of Intralipid 20% (Fresenius Kabi, s-Hertogenbosch, Netherlands). 
The reduced scattering coefficient of the Intralipid stock solution was determined to be 
/ij(800nm) = 18.1 mm by using the formula given by Michels et al. [36]. The depen- 
dence of jj.' s on the intralipid concentration was assumed to be linear and phantoms with 
^(800nm) = [0.072,0.144,0.24,0.36,0.48,0.72,0.96, 1.2,2.4,3.6] mirr 1 were prepared to a 
volume of 20 ml. Reflectance measurements were performed by submerging the probe tip a few 
millimeters into the phantom and suspending it multiple centimeters from the container bottom 
or sides, such that the probe would not collect reflections from container surfaces. The measure- 
ment was repeated 3 times for each phantom, with these spectra averaged for each phantom. 
All measurements were performed with 5 fiber probes varying in diameter: [0.2, 0.4, 0.6, 0.8, 1] 
mm. The phantom was gently agitated before collecting each spectrum to ensure homogeneity. 



2.6. Absolute calibration of experimental single fiber reflectance spectra 

Measurement of single fiber reflectance intensity from the set of solid Ti02 phantoms (outlined 
in Section 2.4) yielded R™p' s ~ rd vs. Ti0 2 weight percentage (wt %) for a range of Ti0 2 wt % 
and fiber diameters. Ti02 wt % was converted to fl' s by estimating a linear conversion factor. In 
order to do this, measurements on the set of solid phantoms from all fibers were interpolated 
onto a curve that was normalized to the highest measured intensity in the set; this conversion 
corrected for different probe-spectralon distances, as noted in Section 2.3. Then a linear con- 
version factor (n' s I Ti02 wt %) was estimated by scaling the normalized reflectance vs. Ti0 2 
wt % data to fit the reflectance predicted by the mathematical model, presented in Eq. (2). This 
calculation yields a jj.' s for each Ti02 phantom. 

Measurements of liquid Intralipid phantoms were converted from a value relative to spec- 
tralon distance, R m p' s ~ rel (IL), to an absolute percentage scale, R^. as ~' 



R' 



meas—abs 
SF 



(IL) = R™p as ~ rel (IL) 



R 



meas—abs 
SF 



s—rel 



l SF 



(phantom) 
(phantom) 



(6) 



where R 



meas—rel 
SF 



(phantom) is the relative reflectance intensity measured from a solid Ti02 phan- 
tom after the calibration procedure given in Eq. (5) and Rf F as -" bsl (phantom) is the absolute 
reflectance on that phantom expressed as a percentage of incident photons, which is calculated 
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from the known /ij of the phantom and the mathematical model of reflectance given in Eq. (2) 
(using values for [pi,p2,p3] specified in Section 3.2.1) . 

3. Results 

3.1. Mathematical characterization ofR^f dependence on scattering properties 
3.1.1. Influence of phase function on R^f vs. dimensionless scattering 

Figure 3 shows the MC simulated reflectance intensity (Rff) vs. dimensionless scattering 
(jj.' s dfib) with combinations of a wide range of fl' s [0.1 — 300] mm -1 and d/n, [0.2 — 1.0] mm. 
Here, Rff is expressed as a percentage of the incident photons. In Fig. 3, symbols differen- 
tiate the specified phase function, including R¥£ data from the HG (g values in the range of 
[0.5 — 0.9]) and MHG (effective g = 0.9) phase functions (see legend). Inspection of the data 
for each investigated phase function, shows that Rff increases monotonically vs. fJ-' s dfn, until 
approaching an asymptotic limit for high dimensionless scattering values. For all phase func- 
tions investigated, the R^f data show a transition from a phase-function dependent regime at 
low n' s d fa, to a region insensitive to phase function at higher }i' s d with the transition point 
between the two regimes occurring for all phase functions at approximately )J.' s dfib > 10. In the 
low dimensionless scattering region, there is a clear stratification between reflectance intensities 
measured from different phase functions, with increased reflectance attributable to an increase 
in the backscattering component of the phase function. Therefore, for the same effective g, Rff 
is greater for the MHG than for the HG phase function; the difference is attributable to more 
high angle backscatter events in the MHG phase function, which is evident from the plot of 
the phase functions displayed in Fig. 2. Moreover, the Rff from HG simulations stratifies in 
order of decreasing backscatter component, with greater reflectance collected from phase func- 
tions containing a higher percentage of large angle backscatter events (e.g. g = 0.5 resulted in 
a higher reflectance than g = 0.9). Conversely, in the region of high dimensionless scattering 
(jl'gdfn, > 10) the Rff shows an insensitivity to backscattering features of the phase function, 
as data from all phase functions converge until they overlap and then asymptotically approach 
a limiting reflectance value. These results indicate that Rff is dependent on the phase function 
in a region of low-mid dimensionless scattering but approaches phase function independent 
behavior in a region of very high values of n' s dfib- 
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Fig. 3. Single fiber reflectance vs dimensionless scattering for HG phase function 
(g=[0.5,0.7,0.8,0.9]) and MHG phase function (g=[0.9]). 
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3.1.2. Mathematical model of R^f vs. dimensionless scattering 

The observed relationship between Rgf and n' s dfn, was described using the mathematical 
model introduced in Section 2.2. Figure 4(A) shows R^f for the HG phase function with g = 0.8 
for /i£ = [0.1 — 300] mm -1 and dfu, = [0.2 — 1.0] mm (blue cross marks), and the correspond- 
ing model predicted reflectance (R^p del ) estimated from Eq. (2) (black line). The estimated 
parameter values, given in Table 1, result in predictions of R^ del that are highly correlated 
with simulated reflectance intensities over the full dimensionless scattering range, as evidenced 
by a Pearson Correlation Coefficient of r = 0.999 and a mean residual error of < 3%. 



(A) (B) 




icr z 10° , io 1 10 2 10 3 



Fig. 4. (A) Single fiber reflectance and mathematical model fit for HG phase function 
(g=0.8). Model estimates of incident photons contacting fiber face (B) and collection effi- 
ciency of fiber (C). 



Equation (2) expresses Rfp as the product of O, the fraction of incident photons that are 
remitted from the medium in contact with the probe face, and 7] c , the fraction of photons that 
contact the fiber face within the cone of acceptance and are collected. Figures 4(B) and 4(C) 
show the MC model estimates (green squares and red diamonds) and mathematical model pre- 
dictions (black lines) of O and t] c , respectively, expressed as percentages, as calculated from 
Eqs. (3) and (4). Figure 4(B) shows that O increases with fl^dm and exhibits a saturating be- 
havior as it approaches the 100% upper limit, a feature well-described by the model. Figure 
4(C) shows that T] c shows a decreasing trend vs. dimensionless scattering, until approximately 
)J-' s dfjb < 1, at which point it approximates a limit linked to the fiber opening angle of accep- 
tance; this behavior is consistent with a previous theoretical analysis [33] and is discussed in 
detail in Section 4. It is important to note that the estimated model parameters in Eq. (2) were 
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not fit to MC estimates of O and r) c ., as they are not observable experimentally. Instead O and r\ c 
were calculated from the parameters resulting from the fit of R^p del from Eq. (2) to the data, as 
shown in Fig. 4(A). The results presented in Figs. 4(B) and 4(C) suggest that the mathematical 
expressions used to represent O and r\ c as components of Rff del are, in fact, representative of 
the true underlying factors contributing to collected reflectance. 

Figure 5 shows R^' del vs. Rff for jj.' s = [0.1 - 300] mm- 1 and d fih = [0.2 - 1.0] mm, 
with symbols distinguishing each phase function, including HG (g = [0.5 —0.9]) and MHG 
(g = [0.9]). Model parameters for Rfp del were fit for each simulated phase function; estimated 
parameter values for each fit are shown in Table 1 (top). Model predictions show an excellent 
agreement with simulated reflectance intensities over the range of phase functions investigated, 
with all r values above 0.998; additionally, the mean residual between all model estimated and 
MC simulated data points was < 3%. These data indicate that the dependence of R^ del on 
dimensionless scattering can be accurately estimated for a given phase function. 




Fig. 5. MC simulated single fiber reflectance vs. model estimate. 



3.2. Experimental determination of jl' s from R$f 
3.2.1. R'^ as -" bs measured from solid Ti0 2 phantoms 

Reflectance measurements on TiCh phantoms yielded R%£, as ~ rel vs. TiC<2 wt % data in the range 
[2.5 - 30] wt %. As outlined in Section 2.6, these data were scaled to R™ as - abs us i n g the 
semi-empirical reflectance model. Because the exact phase function of the TiC»2 phantoms is 
unknown, we used the averages of the fitted parameter values: pi = 7.9, p2 = 1.08, p$ = 2.4. 
Using these parameter values, a linear conversion factor of (jj.' s I TiCh wt %) = 7.42 ±0.35 
was found. Figure 6 shows the R™ as ~ abl '(800nm) from all measurements on TiC>2 phantoms 
vs. jJi' s dfib as measured with 5 probes with diameters in the range [0.2— 1.0] mm, using the 
7.42 (/i'/ TiC<2 wt %) ratio. Also displayed on the plot are Rff vs \l' s dfa> from simulations of 
an HG phase function with g = 0.8, which shows good agreement between experimental and 
simulated data. Note that due to the high reduced scattering coefficients of the phantoms, these 
measurements are in the high dimensionless scattering region (^dfn, > 10), which showed an 
insensitivity to phase function. Varying the parameters pi — p3 in Eq. (2) to be representative of 
HG (g = [0.5 - 0.9]) and MHG (g = [0.9]) phase functions (see Table 1) induced only a 2.2% 
mean residual error between R™ as - abs (based on the average p\ — p^) and Rfp del (based on 
the different p\ — pi combinations listed in Table 1). These results led to the selection of the 
30% Ti02 phantom (which equated to \i' s = 222.6 ± 10.5 mm -1 ) as the calibration phantom for 
subsequent measurements on Intralipid. 
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Fig. 6. Single fiber reflectance measured on TiCh phantoms and simulated by MC model 
(HG phase function with g=0.8). 



3.2.2. R™ as abs measured from Intralipid phantoms 

Figure 7(A) shows R™p is ~" bs from measurements of Intralipid at multiple wavelengths (A = 
[400 — 900] nm in steps of 100 nm) vs. dimensionless scattering. Here the R™p is ~" bs values 
represent the raw intensities scaled by measurements on the solid calibration phantom, allow- 
ing the conversion to an absolute scale. This calibration technique allows reflectance intensities 
from measurements with 5 fibers with different diameters to fall onto the same curve. Figure 
7(B) shows Rff vs. n' s dfu, from MC simulations of the reflectance measurements in Intralipid 
over the same range of wavelengths. Here, the phase function for each wavelength was specified 
using an empirical expression (as in [36]), and the corresponding anisotropy values were in the 
range g = [0.477 — 0.818]. For both measurements and simulations, the R$f values measured 
at different wavelengths (and therefore different phase functions) show stratification at low di- 
mensionless scattering values, and gradually collapse towards a limiting value; these features 
are similar to the data shown in Fig. 3 for HG and MHG phase functions. Figure 8(A) dis- 
plays the corresponding simulated Rff vs. measured R™p ls ~" bs values. These data are highly 
correlated (r = 0.998) and show good agreement (with a mean residual of < 8%) over a range 
of measured reflectance values that span 2 orders of magnitude (range = [0.0158 — 1.57]%). 
It should be noted that relative deviation between simulation and experimental data increases 
for decreasing measured reflectance, and in turn for decreasing dimensionless scattering; this is 
discussed in Section 4. After accounting for the index of refraction mismatch between experi- 
mental and simulated phantoms (1.33 vs. 1.37, respectively), a linear fit of these data through 
the origin gives a slope of 1 .04; potential reasons for the slight offset are given in Section 4. 

3.2.3. Experimental estimation of n' s from I^f u ~ abs 

The semi-empirical reflectance model given in Eq. (2) was fitted to Rff simulated in Intralipid 
for each wavelength-dependent phase function. Table 1 (bottom) lists detailed data for each 
simulation, including the anisotropy, estimated parameters values, and Pearson correlation co- 
efficients for each fit. These model fits were then used to analyze R^, as ~ abs measured exper- 
imentally from Intralipid phantoms and estimate ji' s . Prior to this calculation, the R™ as ~ abs 
data are corrected for slight offset with MC simulations, as noted in Section 3.2.2. Figure 8(B) 
shows /i' estimated using /?^ as ~ and the mathematical model vs. the known experimental 
value; these data are highly correlated (r = 0.997) and show good agreement, with a mean 
residual of 9%. As observed in the reflectance data, the relative error in fl' s estimates increases 
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Fig. 7. Single fiber reflectance measured (A) and simulated (B) from Intralipid optical phan- 
toms at multiple wavelengths X = [400,500,600,800,900] nm; other wavelengths meas- 
ured in this range not shown to improve clarity. 



as dimensionless scattering decreases; this is discussed in Section 4. 

(A) (B) 




Fig. 8. Comparison of (A) single fiber reflectance intensity from MC simulations with ex- 
perimental measurements from Intralipid phantoms A = [400 — 900] nm; (B) jx[ estimated 
from R'Z m ' abs vs. known value. 



4. Discussion 

This paper presents a methodology to quantitatively extract jj.' s from the reflectance intensity 
collected by a single fiber reflectance device (Rsf)- Monte Carlo simulations were used to in- 
vestigate the dependence of Rsf on the scattering properties of an optically sampled turbid 
medium. Simulated data were used to identify a relationship between Rgp and /J.' s dfn„ which 
showed that Rsf approached phase function independent behavior for n' s dfib > 10. This the- 
ory motivated the construction of a solid optical phantom with a very high reduced scattering 
coefficient (fl' s > 200 mm -1 ); such that measurements of Rsf from this phantom are indepen- 
dent of phase function and can be used to calibrate any measured reflectance intensity to an 
absolute scale. Experimental data show that this methodology can be used to accurately extract 
jl' s from Rsf measurements of optical phantoms over a wide range of n' s = [0.1 — 8] mm -1 
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and dfib = [0.2 — 1] mm for a specified phase function. The work presented here provides an 
approach to quantitative estimation of both jj.' s and jj. a from turbid media, such as tissue. 



Table 1 . Estimated Parameter Values Resulting from Fits of Single Fiber Reflectance Model 
(Eqs. (2)-(4)) to R SF Simulated by Monte Carlo Models for HG and MHG Phase Function 
(top) and Wavelength-Dependent Phase Functions [36] in Intralipid (bottom) 



Phase Function 


g 


Pi 


P2 


P3 


r 






Value 


CI 


Value 


CI 


Value 


CI 




HG 


0.9 


9.374 


0.406 


1.152 


0.027 


2.338 


0.851 


0.9994 


HG 


0.8 


8.866 


0.321 


1.120 


0.022 


2.295 


0.667 


0.9995 


HG 


0.7 


8.315 


0.270 


1.086 


0.020 


2.267 


0.561 


0.9996 


HG 


0.5 


7.251 


0.199 


1.018 


0.017 


2.108 


0.405 


0.9998 


MHG 


0.9 


5.676 


0.383 


1.034 


0.046 


3.125 


1.131 


0.9987 


Phase Function 


g 


Pi 


P2 


P3 


r 






Value 


CI 


Value 


CI 


Value 


CI 




Intralipid (400nm) 


0.818 


9.496 


0.103 


1.245 


0.011 


2.995 


0.362 


0.9998 


Intralipid (450nm) 


0.783 


10.001 


0.093 


1.208 


0.010 


3.144 


0.313 


0.9999 


Intralipid (500nm) 


0.749 


9.927 


0.124 


1.170 


0.014 


3.072 


0.401 


0.9997 


Intralipid (550nm) 


0.715 


9.512 


0.124 


1.127 


0.015 


2.745 


0.375 


0.9996 


Intralipid (600nm) 


0.681 


9.014 


0.119 


1.090 


0.017 


2.533 


0.350 


0.9995 


Intralipid (650nm) 


0.647 


8.389 


0.091 


1.054 


0.015 


2.170 


0.255 


0.9996 


Intralipid (700nm) 


0.613 


7.816 


0.077 


1.021 


0.014 


2.010 


0.218 


0.9996 


Intralipid (750nm) 


0.579 


7.361 


0.077 


0.983 


0.016 


1.723 


0.209 


0.9995 


Intralipid (800nm) 


0.545 


6.899 


0.077 


0.947 


0.018 


1.437 


0.204 


0.9994 


Intralipid (850nm) 


0.511 


6.592 


0.060 


0.935 


0.016 


1.461 


0.168 


0.9997 


Intralipid (900nm) 


0.477 


6.323 


0.054 


0.911 


0.016 


1.390 


0.156 


0.9997 



Monte Carlo simulations were used to investigate the effect of jj,' s , dfib, and phase function 
on Rgf . Inspection of the data showed a relationship between R$f and fl' s dfn,, which was char- 
acterized by two regimes: (1) a low dimensionless scattering regime where Rgp depends on 
phase function, and (2) a high dimensionless scattering regime where R$p becomes indepen- 
dent of phase function and asymptotically approaches a limit. The transition point between 
these regimes was observed to be approximately fJ.' s dfji, w 10. The low dimensionless scattering 
regime contains reduced scattering values that are experienced during measurements of tissue 
(common range jj.' s = [0.5 — 5.0] mm -1 ) by single fibers with diameters commonly used clini- 
cally (common range dm = [0.2 — 1.0] mm). While the high dimensionless scattering regime 
will not be experienced in biological tissue, the insensitivity of R$f to phase function in this 
region allows measurements of tissue (of any phase function) to be calibrated onto an absolute 
scale. 

The experimental application of this theory motivated the construction of a solid optical 
phantom with a large concentration of scattering material (TiC»2) that corresponded to jj.' s values 
that were within the high dimensionless scattering (and phase function insensitive) regime. 
These solid phantoms allow single fiber reflectance intensities measured from different fibers 
to be calibrated onto the same absolute scale; previously, this was difficult due to the fiber- 
diameter dependent effect of and variations in probe-spectralon distance utilized in our previous 
calibration technique (as noted in Section 2.3). In this study reflectance measurements on the 
set of phantoms investigated were scaled to a reflectance model to determine the relationship 
between TiC<2 wt % and [i' s . We have assumed that this conversion factor is linear. While it 
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is possible that the relation may not be linear, the data shown on Fig. 6 present a convincing 
argument that this assumption is appropriate. Uncertainty associated with the fl' s / Ti02 wt % 
ratio is not expected to introduce error. This is because the highest scattering phantom (30 wt 
% Ti02 ) was selected for scaling subsequent measurements (on liquid phantoms in this study, 
or on tissue in future work), and in this dimensionless scattering region, a mis-estimation of 
H' s for this phantom are not expected be significant due to the relatively small changes in R$f 
corresponding to changes in jj.' s in this regime. 

The semi-empirical model function used to describe the relation between R$f and dimen- 
sionless scattering for each investigated phase function contains 3 fitted parameters, given by 
Eqs. (2)-(4). The development of this mathematical model provided insight into the factors con- 
tributing to the collected signal. The model represents the light collected during measurement 
as a product of the photons contacting the fiber surface (O) and the fraction of those photons 
within the fiber cone of acceptance (tj c ). Figure 4(B) shows that O increases with dimensionless 
scattering, until it saturates as it approaches the 100% limit. This behavior was described using 
a saturating mathematical expression, which included a phase-function specific fitted power 
(p2) on the dimensionless scattering term; while the fitted values for this parameter were rela- 
tively close to unity (range=[l .018 — 1 . 152], see Table 1) the model fit quality was significantly 
increased by fitting the parameter. Figure 4(C) shows that rj c is large for small dimensionless 
scattering values and decreases as it approaches a limiting value. This transition can be ratio- 
nalized by considering the effect of the optical properties of the turbid medium on the angular 
distribution of remitted light; this concept has been described in detail previously by Bargo et 
al. [33]. In brief, small dimensionless scattering values indicate that photons travel relatively 
large distances between scattering events. Therefore photons remitted from the medium at the 
fiber face have a high probability of being scattered from (relatively) deep locations in the 
medium, and in turn are traveling at angles near normal to the tissue interface. As the dimen- 
sionless scattering increases, photons are more likely to be remitted from shallow distances into 
the medium with remitted light exiting diffusely at all angles. Theoretically, r\ c in this regime 
approaches a diffuse limit defined by the projection of the solid angle of the acceptance cone 
onto the medium/fiber plane divided by all possible solid angles for light remitted across that 
plane; the mathematics underlying this theory has been described in detail by Horn et al. [41], 
and applied to single fiber reflectance spectroscopy by Bargo et al. [33]. It should be noted that 
in the highly scattering region (jj.' s dfib > 10), the angular distribution of remitted light calculated 
by the MC model is not perfectly diffuse; instead, a subcomponent of the light collected by the 
single fiber device has undergone only a few scattering events, which causes slight variations 
(approximately 5%) in r\ c for different phase functions. This phenomenon is known to occur 
for optical devices with overlapping source and detection areas on a medium interface, and 
has been described in detail previously by Snabre et al. [42]. As noted in Section 2.2, Eq. (4) 
is dependent on the fiber NA. Simulations investigating variations of NA on R$p showed that 
the effect could be characterized by adjusting T\u m u exclusively (without refitting parameters 
[pi,p2,p3]), with w 5% mean error between estimates of Rsf measured by fibers of NA= [0.22] 
and NA= [0.1,0.4] in the biologically relevant scattering range (jj.' s = [0.44 — 4]mm~ 1 ) over a 
range of phase functions, with increasing deviations associated with decreasing \i' s d fib (data not 
shown). 

This study includes an experimental proof of concept for the extraction of jJ.' s from experi- 
mental measurements of optical phantoms containing Intralipid. As part of that investigation, 
Rsf was simulated by a MC model that emulated the experimental measurements in optical 
phantoms. These techniques yielded Rsf values that were highly correlated, as observed in 
Fig. 8(A) and 8(B). However, as noted in the Results, there was a slight discrepancy between 
the simulated and the measured Rsf from Intralipid phantoms, with a 4% increase in meas- 
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ured value vs. the corresponding simulated values. One explanation for this slight mismatch is 
that the probes used for the experimental measurements have the fiber tip polished at an angle 
of 15 degrees, altering the fiber face from a circle to an elliptical cross-section, and increas- 
ing the area by approximately 3%; this mismatch is expected to represent the bulk of the 4% 
offset. Additionally, it was observed that the relative error in R$f estimates increased with de- 
creasing dimensionless scattering. This is most likely caused by slight mismatches between the 
actual phase function of the liquid phantoms and the model approximation utilized in the MC 
model [36]. Another source of mismatch may be in the specification of the index of refraction 
of both fiber and sample in the MC simulation, which was assumed wavelength to be indepen- 
dent. Slight mismatches between experiment and simulation would be observed especially in 
regions of small dimensionless scattering values where the reflectance signal is most sensitive 
to phase function. 

The goal for performing the absolute R$f calibration was to extract an absolute value for 
within the sampled turbid medium from the measured single fiber reflectance intensity. This 
was achieved in this study by performing measurements on Intralipid phantoms, calibrating the 
measured reflectance with measurements on the highly scattering solid calibration phantom, 
and analyzing the data using a phase-function specific mathematical model describing the de- 
pendence of reflectance on dimensionless scattering. The work in this study calculated }i' s for 
5 fiber diameters, on Intralipid phantoms of 10 different concentrations at 11 different wave- 
lengths (range [400 — 900] nm), representing a wide range of phase functions and jj.' s = [0.05 — 8] 
mm . Because this calculation utilized the model fit to MC simulations to analyze measured 
reflectance data, the measurements were corrected for the observed 4% offset prior to es- 
timation, as described in detail earlier in this Section. After this correction we were able to 
accurately extract the reduced scattering coefficient from Intralipid measurements, with a mean 
residual of 9% over the entire range of this fl' s range is wider than experienced in tissue, and 
therefore error associated with biologically relevant fx' s values is expected to be < 9%. In this 
study the semi-empirical model was used in conjunction with measured reflectance on phan- 
toms to provide absolute estimates of jj.' s , estimates that have been demonstrated to be accurate 
for a given phase function. In tissue, the phase function is generally unknown and assumptions 
about the values for pi — p3 that feed into the model expression must be made. Future studies 
will investigate how the assumed parameter values affect the estimates of /i', and whether the 
parameter values can be expressed in terms of metrics extracted from a known wavelength- 
dependent phase function profile. Furthermore, we are currently investigating methods to mea- 
sure the parameter values experimentally for biological tissues. 
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